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A Bianchi Type IV viscous fluid model is analyzed for asymptotic isotropy behaviour. In particu- 
lar, we consider the case of a viscous fluid without heat conduction in conjunction with a non-tilted 
Bianchi cosmology in an orthonormal frame using expansion-normalized variables. The Einstein 
field equations in this approach reduce to an autonomous system of first-order ordinary differential 
equations which because of their high degree of non-linearity are solved using numerical Runge- 
■ Kutta methods. It is shown that for cases where the expansion-normalized bulk viscosity coefficient 

£o dominates over the expansion-normalized shear viscosity coefficient, r/o, the model isotropizes. 
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I. INTRODUCTION 



Spatially homogeneous and anisotropic models of the universe have undergone great study and continue to remain 
amongst the most popular areas of research in cosmology. Early-universe cosmological models for the most part 
have assumed the universe to be spatially homogeneous and anisotropic, with the important exception being the case 
of inhomogenous models such as the LeMaitre-Tolman-Bondi and "Swiss-cheese" models. However, if one begins 
with the idea of the early universe being spatially homogeneous and anistropic, then to transition to the present-day 
Friedmann-LcMaitrc- Robertson- Walker models requires the anisotropy in the former models to decay. The process by 
which this anisotropic decay occurs is arguably the most fundamental property of any early-universe model which aims 
to transition to the present-day models. For example, Belinskii, Khalatnikov, and Lifshitz [l| studied the oscillatory 
approach to a singular point in relativistic cosmologies. Misner Q studied the anistropic decay of the vacuum Bianchi 
Type IX/Mixmaster models of the universe. A very general approach to describing the isotropization of Bianchi 
models was described by Salucci and Fabbri Q . Coley and van den Hoogen [3| studied causal anistropic viscous fluid 
models and described conditions for such models to isotropize. As for very recent work on this subject, Pradhan, Rai, 
and Singh [j| studied the Bianchi Type V bulk viscous models and showed that such models do isotropize for specific 
functional forms of the anistotropic scale factors. 

As for Bianchi Type IV models specifically, Hervik, van den Hoogen, and Coley @ studied future asymptotic 
1 behaviour of tilted vacuum Bianchi Type IV models, and found that such models do not necessarily isotropize at 
. late times. Uggla and Rosquist Q studied the orthogonal Bianchi Type IV model near the initial singularity with 
r* 1 a vacuum or perfect-fluid source. We chose to study the isotropization behaviour of the Bianchi Type IV viscous 
model largely because such a study has not been taken on extensively in the literature, and perhaps, such a study 
will add to the already rich landscape of spatially homogeneous and anisotropic models of the early universe. We 
will employ the important technique of orthonormal frames pioneered by Ellis and MacCallum [8| which reduces the 
Einstein field equations which arc a coupled set of ten hyperbolic nonlinear partial differential equations into a system 
of autonomous nonlinear ordinary differential equations. 

II. THE ORTHONORMAL FRAME APPROACH TO THE BIANCHI TYPE IV ALGEBRA 

If one considers a fluid described by a family of time- like curves, with 4- velocity, u°, the properties of the fluid flow 
are most appropriately described by the decomposition [9| 

Ua;b = -a a U b + <7 ab + LO ab + -h ab 8. (1) 



'Electronic address: isk@yorku.ca 



^Electronic address: mchaslam@mathstat.yorku.ca 



2 



In this equation 9 = u? a is a measure of the divergence of the family of time-like curves, a a = u a: i,u b is a measure 
of how much the curves represent non-geodesics (and thus, can be taken to be the fluid acceleration), and h a {, = 
g a b — u a Ub/u c u c is a projection tensor. We also have terms that represent the shear and vorticity of the fluid: 



0~ab 



1 1 



hTh 



b ' 



(2) 



^ab — ^ "njm) ■ 

In addition, we note that since for any velocity field u a u a = ±1, the following relationships also hold: 



a c u 



(3) 



(4) 



When using orthonormal frames, the metric tensor, g uv takes a very simple form for the following reason: any 
cosmological model can be described by a coordinate system in which we can define a basis {e M }, and its dual basis 
of differential one-forms {u u }- As usual, the metric tensor is then defined as g uv = g(e u ,e v ), from which the line 
element is ds 2 = g uv u> u uj v . For a given coordinate chart on the pseudo-Riemannian manifold, we can take the basis 
formed by {e u } as the coordinate basis {d/dx 1 }, with the dual basis then playing the role of the coordinate differential 
one-forms, {dx 1 }. If we now assume that the basis constitutes an orthonormal frame o, then all four basis vectors 
{e u } are mutually orthonormal, and thus satisfy the relationship: g(e u , e v ) = rj uv = S uv = diag(—l, 1, 1, 1). Hence, 
g uv is really just the Kronecker tensor, and all tensors in this approach have their indices lowered(raised) with it. 
Furthermore, the group structure constants are functions, which we denote by a lowercase c to distinguish between 
the constants typically described in the non-orthonormal frame approach by an uppercase C (Io| (llj . These functions 
satisfy the commutation relation [l2j [H| 



Recall that for an arbitrary basis, for a general (torsion-free) spacetime, wc have 

a pa pa 

uv vu uv ' 

We can then write the connection coefficients as 



(5) 
(6) 



2 (9ab c vu "t" 9uvC av 9vb c au) 



(7) 



Because we arc assuming a non-tilted cosmology, the fluid four-velocity is necessarily orthogonal to the spatial 
hypcrsurfaces, which implies that the fluid is irrotational and geodesic. That is, 



0. 



We can therefore write that 



Exploiting the non-tilted condition once more, we see that Eq. ([9]) implies 



T £ =S> c 

UV 7 



0. 



(8) 
(9) 
(10) 



Using the Behr decomposition, the structure coefficients(functions) can be decomposed into symmetric and anti- 
symmetric parts as 



Hb 



(11) 



The vector O c has the significance of being interpreted as a local angular velocity in the rest frame of an observer 
with spatial frame {e a }. It is conventionally defined as [l2| 



n a = -e abcd 



u b e c ■ e d . 



(12) 
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It can be shown that the structure coefficients are spatial, and therefore, must correspond to the Bianchi models. We 
can therefore write that 



eijii 



Jk 



ai (5lS 



l^k 
3 



(13) 



Since each structure coefficient is constant along each orbit of transitivity, n ab and dj are functions of time only. We 
can find evolution equations for n ab and a, by first realizing that the Jacobi identity holds for all vectors. The Jacobi 
identity as applied to the set of vectors (U,e a ,eb) gives 



[U, [e a , e 6 ]] + [e a , [e b , U]] + [e 6 , [U, ej] = 
* (U(<&) + «6 + c: u cl + cLc^K = 0. 



(14) 



In this derivation, U is a gauge vector that will typically be set to cV Note that since we are assuming that the fluid 
is non-tilted we necessarily have 9 UV = F^, => c\ a = c} ab = 0. We therefore obtain the identity 



dt{c k ab 



+ c td c ab + c . 



adHt 



+ c bd c ta 



0. 



(15) 



One can show that upon applying Jacobi's identity (Eq. (fT4"]) ) to the three spatial vectors, we get the eigenvalue 
equation [l2j 



n lJ a,i 



0. 



The trace of Eq. (|15[) gives the evolution equation for a, as 



-Oai 
3 



(7i 3 o? + e^a-T^ = 0. 



On the other hand, the trace- free part of Eq. ([T5|) gives the evolution equation for n ab as 

1 



n ab 



-9n 



ab 



2n k{a a b ) k = 0. 



(16) 



(17) 



(18) 



The structure constants as defined in Eq. (|13p must correspond to a Lie algebra as per the Jacobi identity (Eq. 
(fl"4|) ). thus, Eq. (fT6|) must hold for any choice of riy and a^. Bianchi models are typically classified into two categories: 



Class A models and Class B models. Class A models are all Bianchi models such that a, = for which Eq. (fTBj) is 
satisfied. For class B models, ai must be an eigenvector of n lJ with zero eigenvalue. It is important to note that we 
always take n lJ to be a symmetric matrix and as such we can diagonalize it using an orientation of our choice for the 
spatial frame. It has been conventional to assume 



n.ij = diag (711,112,113), , a 1 = (0, 0, a). 



(19) 



The Jacobi identity (Eq. (fT4|) ) immediately implies that 71330 = 77,3a = 0. It is a fact of linear algebra that eigenvalues 
of a matrix are invariant under a conjugation operation with respect to rotations. One can then classify the different 
Bianchi models by the signs of the eigenvalues Tin, n%2, ^33 = tii, 712,71-3 and a. These have been listed in Table 1. 
Since, we are interested in the Bianchi Type IV algebra, we choose 
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TABLE I: Bianchi Type classifications in terms of the time- dependent eigenvalues. 



aS l 3 > 0, and ri\ > 0, n 2 = n 3 = 0. 



(20) 



4 



III. THE ENERGY-MOMENTUM TENSOR FOR A VISCOUS FLUID 



In this section, we will derive the form of the energy-momentum tensor under concern, namely, for that of a viscous 
fluid without heat conduction. We will then derive the dynamical equation for this fluid, and will see in the next 
section that the latter closes the system of dynamical equations. Recall that the energy-momentum tensor for a 
perfect fluid takes the form 

T ab = (ji+ p)u a u b - u c u c g ab p. (21) 
For the moment, letting fx + p = w, we obtain 

T ab = wu a u b - u c u c g ab p, (22) 
Denoting the viscous contributions by V a b, we seek a modification of Eq. (f22j) such that 

T ab = wu a u b - u c u c g ab p + Vab- (23) 

To obtain the form of this additional tensor term, we note that from classical fluid mechanics, the Euler equation is 
given as 

{p Ui ) t = -Ihk,k, (24) 
where Iljfc is the momentum flux tensor. Also, recall that for a non- viscous fluid, one has the fundamental relationship 

H-ik = pS lk + puiUk- (25) 
We simply add a term to Eq. (|25[) that represents the viscous momentum flux, S^ fe , to obtain 

Ihk = pSik + pUiUk - = -£ lfe + puiu k . (26) 

It is important to note that 

= -pSik + 14 (27) 

is the stress tensor, while, T,' ik is the viscous stress tensor. The general form of the viscous stress tensor can be formed 
by recalling that viscosity is formed when the fluid particles move with respect to each other at different velocities, so 
this stress tensor can only depend on spatial components of the fluid velocity. We assume that these gradients in the 
velocity are small, so that the momentum tensor only depends on the first derivatives of the velocity in some Taylor 
series expansion. Therefore, Y/ ik is some function of the In addition, when the fluid is in rotation, no internal 
motions of particles can be occurring, so we consider linear combinations of u^k + Uk,ii which clearly vanish for a fluid 
in rotation with some angular velocity, fli. The most general viscous tensor that can be formed is given by 

= V (ui^ k + u k ^ - ^S lk ui^ + £6ikUij, (28) 

where r\ and £ are the coefficients of shear and bulk/second viscosity, respectively [ljj [lj|. In Eq. (|28|) . we note that 
SikUi^i is an expansion rate tensor, and [u^k + "fc,? — § $ikUi,i) represents the shear rate tensor. Since we would like 
to generalize this expression to the general relativistic case, we replace the partial derivatives above with covariant 
derivatives, and the Kroenecker tensor with a more general metric tensor, that is, Sik — > g%k- We thus have that 

Sifc = V ^iii-k + u k . t - ^g ik ui-i^ + £,g ik ui-i. (29) 

Denoting the shear rate tensor as a ab , and the expansion rate scalar as 8 = u a a , Eq. (|29p becomes 

Vab = -2 Wab - ZOhab- (30) 

Substituting Eq. (f3T>l) into Eq. (f2"3"f we finally obtain the required form of the energy-momentum tensor as 



Tab = (/x + P) u a u b - u c u c g a bP - 2Wab - £,9h a b- 



(31) 
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A. An Equation of Motion for the Fluid 

Now that we have the energy- momentum tensor in hand (Eq. pip ), we are in a position to derive an equation of 
motion for the fluid. Equations of motion are typically found from conservation laws. The generalization of conser- 
vation of energy-momentum in general relativity appears via the divergence-free property of the energy-momentum 
tensor, which itself, is coupled to the divergence-free property of the Einstein tensor. For notational convenience 
(which will be clearer in what follows), we will let n a f, = 2?7cr ttb , such that 

T b = ( M + p)u a u b + S b p -n b a -t9 (S b + u a u b ) . (32) 

Energy conservation requires that this quantity should have a vanishing divergence, T.^ = = T b . b . Taking the 
divergence of Eq. (|32[) and contracting the resulting equation with u a , we obtain 

u a (Gu + p)u a u b + S b aP -TT b a -Z6 {6 b a + u a u b )). b = 0. (33) 

We will proceed by evaluating Eq. (|33[) term-by-term. For the first term, we obtain 

u a (([i + p)- b u a u b + (p + p)u a - b u b + (/i + p)u a u\) 



(-fi-p-(v+ P )9). (34) 



For the second term, we obtain 



For the third term, wc have 



N 

p. (35) 



v ab K ab - (36) 



For the fourth term, we have 



u a {S b +u a u b )). b 



= -£6 2 . (37) 

Substituting Eqs. (J34J) - (37]) into Eq. <j33j), we obtain 

A + (a* + P)9 - 2 Wab a ab - iO 2 = 0. (38) 

Using the definition a 1 = \a ab a ab =>■ cr ab u a b = 2cr 2 , we have finally that 

(i+(jj,+p)6- 4rycr 2 - £6 2 = 0. (39) 

Before wc conclude this section, we note the necessity of defining a barotropic equation of state. We will choose 
p = ^fi, which comes from the idea of the fluid moving ultra-relativistically as would be the case in the early universe. 
Therefore, substituting this equation of state into Eq. (j3"9"]l . we obtain 



A + (jA'J o - 4 W 2 - te 2 = 0. (40) 

This is the evolution equation for the fluid as it propagates through the universe. As we shall see in the next section, 
this is the dynamical equation which closes the Einstein field equations. 



IV. THE DYNAMICAL EQUATIONS 



Using the orthonormal frame approach we will now write the Einstein field equations that correspond to the Bianchi 
Type IV algebra. The resulting equations are the Raychaudhuri equation, the generalized Fricdmann equation, and 
shear propagation equations. The system will be closed by Eqs. (JTTJ) , (TT5)) and PO)) . 
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A. The Raychaudhuri Equation 

We now derive the viscous analogue of the Raychaudhuri equation which describes the expansion/contraction 
behaviour of a congruence of fluid lines, and necessarily the universe. We begin by recalling that one can consider a 
fluid evolution as described by a family of time- like curves, with 4- velocity, u a . We had from Eq. ([1} that 

u a -b = ~a a u b + a ab + uj ab + -h ab 9. (41) 

Analyzing the equation of the divergence of the fluid flow, we are interested in how this divergence evolves in time. 
That is, we are interested in the quantity 

e=(u« a ). b u b = u^ b u b . (42) 
By the definition of the Ricmann curvature tensor, we have that 

<a;6 - U% a = -R bc U C . (43) 

Substituting Eq. (J43J) into Eq. (J42J) , we obtain 

= u%. a u b - Rbcu c u b . (44) 
Using the identity {u a b u b ^j = u a b . a u b + u a b u b a , Eq. (|4"4")) becomes 

9 = a% - R bc u c u b - u a;b u b - a . (45) 
Since h ab h ab = 3,0-2 = \a ab a ab ,u? = \u ab u a b, h ab a ab = 0, and 

h ab u a , b = -a a u b h ab + h ab a ab + h ab uj ab + ^9h ab h ab , (46) 



one obtains 



9 = -R bc u b u c + (a a a - T a an a n ) + 2lo 2 -2<t 2 -^9 2 . (47) 



Rab = k ( T ab - ]-Tg ab I . 



The Einstein field equations imply that 

R .u = k I TLl — 

2 

Computing the trace of the energy- momentum tensor (Eq. (|32|) ) we obtain 

T = u a u a {n - 3p) - 2rjg ab a ab - tg ab (9h ab ) . (49) 
Contracting Eq. © with the metric tensor g ab we get that 



f b <Jab = g ab \ U( a . b ) + a (aUb) - ^9h ab ^j = 0. 



(50) 



In addition, note that the contraction of 9h ab with the metric tensor g ab simply gives 9g ab h ab = 39. Therefore, Eq. 
(|15| now takes the simple form 

T = u a u a (fi-3p) -3£6. (51) 
We also make use of the additional property, namely that 

T ab u a u b = p. (52) 

Substituting Eqs. (gSJ) - (J52J into Eq. (gSJ) we get 

R ab u a u b = i/s (ji + 3p) + ^9u a u a . (53) 
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We will use the metric signature (— 1, +1, +1, +1) so that u a u a = — 1, and Eq. ([53]) becomes 

R ab u a u b = X -k (fi + 3p) - ^9. (54) 

Substituting Eq. (|54[) into Eq. ()47|) we obtain the Raychaudhuri equation for our viscous fluid as 

9 = ~k (fx + 3p) + + < + 2w 2 - 2<r 2 - itf 2 . (55) 

Since we are assuming a non-tilted cosmology, the fluid is necessarily geodesic and irrotational, so that Eq. (|55|) upon 
applying the ultra relativistic equation of state becomes 

9 = -kh+ -k£,9- 2cr 2 - -9 2 . (56) 
2 3 

B. The Generalized Friedmann Equation 

The generalized Friedmann equation is an extension of the standard Friedmann equation one obtains from the FLRW 
metric. This equation is important because it relates the rate of expansion to the shear, and more importantly, the 
curvature of the spatial slices in the spacetime foliation. Exploiting the assumption of a non-tilted cosmology again, 
we state the additional fact that in this case 

Ua:b = 9ab = Kab- (57) 

Here, we have denoted the extrinsic curvature of the spatial slices by the extrinsic curvature tensor, K ab . Recall from 
standard differential geometry, we have that for any three-dimensional spatial slice, 

= ^R + K 2 - K a ?K a0 - 2^R a pu a uP 

^ K T a pu a u = \( {3) R-K a0 K aP + K 2 ). (58) 



2 

The decomposition equation (Eq. (Q3) then becomes 



3 a ji = K a p = a af} + -h al3 9. (59) 



Substituting the extrinsic curvature scalar as K 2 = ^K a fjK a/3 , and Eq. (|52p into Eq. ([55| . one obtains the generalized 
Friedmann equation, 

11 1 (3) 

3# 2 = 7fabO- ah - g R + KfX. (60) 

C. The Shear Propagation Equations 

The last set of dynamical equations implied by The Einstein Field equations are the shear propagation equations. 
They essentially describe the evolution of the anisotropy in a cosmological model as a function of time. The lengthy 
derivation of these equations has been done in much of the literature on cosmological dynamical systems [l3[ [IH 
[l6[ ■ The idea behind the derivation is similar to that of the derivation of the Raychaudhuri equation. The shear 
propagation equations are 

o-ab + 9a ab - a d a e bcd n c - a<te acd n c +< 3 ) R ab - ±6$R = 2 V Ka ab . (61) 
D. The Ricci Curvature and Constraint Equations 



Perhaps the greatest accomplishment of the orthonormal frame approach is that the curvature tensors, namely, 
the Ricci tensor are no longer expressed in terms of the coordinate basis functions which simplifies computations 
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considerably. The Ricci tensor and scalar are expressed entirely in terms of a 1 and n a b- As discussed by Ellis and 
MacCallum Q and Gr0n and Hervik [l2j , the Ricci tensor takes the form 

{3) R a b = -e c a d n bc a d - e c b d n ac a d + 

( 1 \ (62) 

2n ad n d - nn ab - 5 a b f 2a 2 + n cd n cd - -n 2 J . 

Contracting Eq. (|62p . we obtain the Ricci scalar, 

MR Rl = ^ 6 a 2 + n cd n cd - l -n^j . (63) 

One can also show that the off-diagonal components of the four-dimensional Ricci tensor yield a constraint equation, 

3a b a ba - e abc n cd a b d = 0. (64) 

For computational purposes, it is more convenient to express Eq. (|64p in component form, 

3aa 33 + (n 11 - n 22 ) a 21 = 
3aa 3 i + n 22 a 3 2 = 

3ao- 32 - n 11 a 3 i = 0. (65) 

Up to this point we have derived and discussed all of the necessary formalism required to adequately describe the 
dynamics of a spatially homogeneous and anisotropic spacetime. In the next section, we will apply these results to 
the Bianchi IV algebra. 

V. THE BIANCHI TYPE IV DYNAMICAL EQUATIONS 

In this section, we will derive the dynamical equations associated with the Bianchi Type IV algebra, and study 
their solution. For the Bianchi Type IV geometry, recall from Eq. (f20|) . we take n\\ = N, and 1222 = "33 = 0, where 
N > 0. The constraint relations Eq. (f65|) then imply that 

3acr 33 + N&2\ = 
3a(73i = 

3acr 32 - 031 = 0. (66) 

We can see from Eq. (f6"6")l that 032 = 0. Since the shear tensor is taken to be symmetric and traceless, in three 
dimensions, it only has two independent components. It is also clear from Eq. (|66[) . that the only non-zero off- 
diagonal component of the shear tensor is a 21 = <J\2- Therefore, for the case where a = b, 



\/3cr_,cr+ - \/3<7_,-2cr + ) . (67) 



The non-zero, independent, off-diagonal component for the shear tensor is given by: 

-3oa- 33 6aa + 

a 21 =a 12 = ^— = —. (68) 

Substituting Eqs. ([67]). ([55]1. ([H]) into Eqs. (JT7J) and ([T5|l. we obtain the relations, 

n 1 =n 2 = 0, O 3 = -6acr+, (69) 

a + \e a - 2a+a = 0, (70) 



and 

N - 

3 



-leN -2N + Vicr-^j =0. (71) 



Using Eqs. (|62|) and ((20)) , we compute the three-dimensional Ricci tensor as 

{3) R ab = -efn bl {a) - e\ d n al {a) + 2n al n\ - Nn ab - S ab (2a 2 + ^N 2 ^j 
Eq. ([72} then gives the components of the Ricci tensor as: 



(72) 



(3) i?n = iiV 2 -2a 2 , (73) 

^R 12 ^ R 21 = Na, (74) 

WR 13 =< 3 > i? 3 i = 0, (75) 

(3) i? 22 = -(2a 2 + liV 2 ), (76) 

^R 23 = (3) i? 32 = 0, (77) 

(3) i? 33 = -(20? + ^. (78) 
The Ricci scalar is computed from Eqs. ([B"3"} and ([20"} as 

i?= (3) R a a = -6a 2 - ^N 2 . (79) 

From Eqs. ([61} . (|73l) - ([78|) . and ([79} . we see that the diagonal components of the shear propagation equations are: 

/ \ 2iV 2 / \ 

cr' + + V3V_ + 9 + V3(T_J - 72a 2 cr 2 _ + = 2r)K (ct+ + x^<7-) , (80) 

/ \ A^ 2 / 

a' + — V3(f- + 9 \ a+ — \/3cr_ J + 72a 2 er 2 _ — = 2r)K [a + — \fStj- 



and 



N 2 

— = 2r)Ka+. (82) 



G 

The off-diagonal component of the shear propagation equations is evidently 

6aa + N r 2 

6der+ + 6a<j + — h Q9aa + + 12v 3aer+c7_ + N a = 12rjnaa + . 



(83) 



Substituting Eqs. (|7Tj) and ([70]) into Eq. ([53"} . we obtain upon simplification 

iV 2 + 6cr'+ + 2cr + 6» = 12«;r/cr + . (84) 
Subtracting Eq. ([50} from Eq. (|5T} . we obtain 

iV 4 - 144a 2 cr 2 _ + 2VZN 2 [a- + <t_ (9 - 2kt))} = 0. (85) 
Adding Eqs. ((52} and ([51}. we obtain 

-5A^ 2 
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5cr'+ - a+9 = -IO/CT70+. (86) 



One can see that Eqs. ([85} and ([86} constitute the dynamical equations for the shear variables <r+ and er_. 
Substituting Eq. ([79"} into Eq. ([50} we obtain the appropriate form of the generalized Friedmann equation as 



a 2 = 1-9 2 - 3a 2 - \n 2 - kll. (87) 
3 4 v ' 



We can use Eq. ([57} to rewrite the dynamical equation for the fluid density (Eq. (00}) as 

fi = -I2a 2 v - An m - V N 2 - ^ + ^f- + 9 2 ^ 
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and the Raychaudhuri equation as 



9 = - (12a 2 + 2k(x + N 2 - 29 2 + 3k(9£) . 

Therefore, one can now see that our Bianchi IV viscous fluid model dynamics is governed by Eqs. (I89[) . 
, dZOD, and (EH): 



(89) 



13, 



- (12a 2 + 2n[i + N 2 - 26 2 + 3n0£) 



-\2a 2 rj — AKfiri — rjN 
N 2 24 v / 3a 2 a 2 



2 4fi6 Ar,9 2 
3 3 



a- = — 



2VS iV 2 
N 2 a + 6 



— (T-9 + 2k(t_?7 



6 5 
-af--2a + 



2kt]g + 



N = N 



-+2[o + + V3<T- 



(90) 



This is a coupled first-order system of six nonlinear ordinary differential equations with free parameters being the 
shear and bulk viscosity coefficients. Clearly, the system (Eqs. (|5U)) ) has no exact solution, so numerical methods must 
be applied. However, because of the high degree of nonlinearity of the system, numerical algorithms can be difficult 
to employ to obtain any relevant solutions. Therefore, we will write these equations in the expansion-normalized form 
fl3j [l2| which will reduce the number of equations to five. 



VI. DYNAMICAL EQUATIONS IN EXPANSION-NORMALIZED VARIABLES 



The basic idea of the expansion- normalized variables is that the system (Eqs. (19010 is of the form 

dx. . , , , 

Tt = (91) 

where x = (9, fj,, er_, <r+, a, N) £ R 6 . Wc will define 9 = 3H, where H is denoted the Hubble scalar and is defined as 

H=-, (92) 
s 

and s is a cosmological length-scale function. It is also convenient to define the cosmological deceleration parameter q 

q=~. (93) 

Clearly, we have the relationship 

H = -(l + q)H 2 . (94) 
It is also necessary to introduce a dimensionless time variable r as 

s = s e T . (95) 

From these equations, one can show that 

Substituting 9 = 3H into the Raychaudhuri equation (Eq. ([89])). wc obtain 

H= -l Kfi + KtH- ^<j 2 -H 2 . (97) 
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We will in addition define: a density parameter: = 3^2, a scalar shear parameter: E 2 = ^jyy, a curvature parameter: 
K = — 4tt§, curvature parameters:?! = N/H, A = a/H, and shear parameters: E± = a±/H. 



Wh 2 

Comparing Eq. (|M|) and Eq. (j9"T)) . we immediately see that 



g = K f]_^_2S 2 . (98) 
H 

At this point we assume that the bulk viscosity coefficient obeys the equation of state 

I = 3£o, (99) 
where £0 is a nonnegative real number, and, Eq. ([98jl now becomes 

q = ntl - 3< - 2E 2 . (100) 

We further proceed by writing Eqs. (|85[) . (|86[) . (|T0|) . and (|71l) in expansion- normalized form using the aforementioned 
definitions of the shear and curvature parameters. Beginning with Eq. (|85p . we see that 

n 2 H n A 2 Y?,H H , v 

£_ = =+2^ ± 3£_iJ + 2«;E_77 £_. (101 

2\/3 n 2 # 

Using the chain rule, we obtain from Eq. (|M|) that 

c2E_ dE_ di ■ 1 , x 

Substituting Eq. ([T02]) into Eq. (fToTj) . we get that 

n 2 A 2 Y? 

£' = = + 24V3 — ^± -3E_ - 6kS_?7o + (1 + <?)E_. (103) 

2^ « 2 

As is done by Saha and Rikhvitsky, [l7| . we have set in Eq. (|103p . T]/H = 3?7o, where 770 is a nonnegative real number. 
Continuing in a similar way with Eq. (|86[) . we have that 

= -In 2 - |e+ + 6kjk,E+ + (1 + g)£+. (104) 
6 5 

In addition, the expansion-normalized forms of Eqs. (|71j) and (|70|) become 

n' = n + 2n ^E + + V3E_) +(l + q)n, (105) 

and 

A' = -A + 2E+V1 + (1 + q) A. (106) 



Note that we calculate E 2 = as 



^ 2 = k772 ( cr afc cra6 ) 



6# 2 
1 



6ff 



2 "T" "22 -T "33; 



V 2 , + 3cr 2 + a\ + 3cr 2 + 4a 2 .) 



6H 

^2 1 v^2 



£i_ + £i. (107) 



In addition, K = — 4tjt is computed to be 



Mr 



K = - 



m 2 



J_ ( Qa 2 + -N 2 
6H 2 \ 2 

^ 2 + ^ (108) 
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The ft term in the expansion- normalized equations is governed by the expansion- normalized version of Eq. (|88l) . 
which turns out to be 

ft' = —12A 2 rjo - UntlTio - ijon 2 - 4ft + 12?7o + 9£ + 2ft(l + q). (109) 

It is important to note that on physical grounds, ft > 0, which acts as a constraint on the dynamical system. This is 
because we assume that there is no exotic matter in our universe model, that is, the matter content is always assumed 
to have a nonnegative energy density. 

In summary, the expansion- normalized dynamical system is governed by Eqs. (|103j) . (|104j) . ()105p . (|106[) . and (|109|) : 

n 2 A 2 Yj 2 

E'_ = = +24^3 — - 3E_ - 6kE_?7o + (1 + g)E_ 

2V3 n 2 

E' + = -in 2 -^E + + 6^ S+ + (l + 9)S+ 
6 5 

A' = -A + 2T 1+ A+{\ + q)A 
n = n + 2n I E + 

ft' = % (-12A 2 - I2nfl - n 2 ) - AVL + 12r/ + 9^ + 20(1 + q). 

(110) 



VII. SOLUTIONS TO THE BIANCHI TYPE IV DYNAMICAL EQUATIONS 

In this chapter, we present several numerical solutions to the system of equations Eq. (jllOp . Clearly, this system, 
which is nonlinear and coupled has no exact solutions. Our goal however is to discover values for £o and 770 such that 
the system asymptotically isotropizes, that is, E± — > as r — > 0. Such physical characteristics resemble an anistropic 
early-universe that over time becomes isotropic. The solutions to this system of equations were obtained by using the 
ODE23 solver in MATLAB. 

Note that in the diagrams below, asterisks indicate the initial conditions for the numerical experiments. We 
simulated solutions to the system of equations Eq. (|110p . by considering ten sets of initial conditions, which had the 
form 

[E_(0),E + (0),A(0),n(0),ft(0)] = [0.05i, 0.07i, 0.02i, 0.03i, O.Obi] , (111) 

with i being enumerated from one to ten. Note that there are many possible choices for initial conditions, however, 
some care must be taken with respect to constraints. In particular, because of the Bianchi IV algebra, one must 
assure that A > 0, n > 0, and Q, > 0. Clearly the initial conditions chosen satisfy these constraints for all values of i 
from one through ten. 
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A. Co = 0.1, 770 = 0.1 



FIG. 1: The figures below display the dynamical system behaviour for £o = 0.1, rjo = 0.1 in terms of the phase plot of the 
anisotropy in the Hubble flow, and plots of the energy density and anisotropy variables as functions of time. One can see that 
fi — > 0, £_ — > 0, but £+ -ft 0, so the model only partially isotropizes. 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



14 



B. 



Co = 0.1, 770 = 1 



FIG. 2: The figures below display the dynamical system behaviour for Co = 0.1, 770 = 1 in terms of the phase plot of the 
anisotropy in the Hubble flow, and plots of the energy density and anisotropy variables as functions of time. One can see that 
O — s- 0, £_ — ¥ 0, but £+ -/-¥ 0, so the model only partially isotropizes. 
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Co = 1, 110 = 0.1 



FIG. 3: The figures below display the dynamical system behaviour for £o = 1, 770 = 0.1 in terms of the phase plot of the 
anisotropy in the Hubble flow, and plots of the energy density and anisotropy variables as functions of time. One can see that 
O — > 0, £_ — > and E+ — > 0, so the model does in fact isotropize. 
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D. 



£o = 10, 770 = 0.1 



FIG. 4: The figures below display the dynamical system behaviour for £0 = 10, 770 = 0.1 in terms of the phase plot of the 
anisotropy in the Hubble flow, and plots of the energy density and anisotropy variables as functions of time. One can see that 
fi — > 0, £_ — > and E+ — > 0, so the model does in fact isotropize. (Note that the time scales of the models were adjusted in 
the figure to display the dynamical behaviour clearly.) 
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E. Limiting Cases of The Viscosity Coefficients 



The numerical analysis thus far naturally leads to the question as to under what conditions does this cosmological 
model isotropize. In an attempt to answer this question, we will take two limiting cases. The first case will be 
for which the expansion-normalized bulk viscosity coefficient is very small, that is, £o ~ and we will let the 
expansion-normalized shear viscosity coefficient, 770 vary. It is therefore appropriate to say in this case that the 
expansion- normalized shear-viscosity is the dominant viscosity coefficient. 



FIG. 5: The different plots in this figure demonstrate the behaviour of the anisotropy functions and energy density as the 
expansion-normalized shear viscosity coefficient is increased while the expansion-normalized bulk- viscosity coefficient is assumed 
to be negligible. One sees that although in each case the energy density is bounded such that Q > 0, E± -fr 0, so the model 
does not isotropize. In fact, one can see that the anisotropic variables "freeze in" to non-zero asymptotic values after a period 
of time. 
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The second case that we consider will be for a negligible expansion-normalized shear viscosity coefficient, that is, 
?7o ~ while the expansion-normalized bulk-viscosity coefficient, £ is allowed to vary. It is therefore appropriate to 
say in this case that the expansion- normalized shear- viscosity is the dominant viscosity coefficient. 



FIG. 6: The different plots in this figure demonstrate the behaviour of the anisotropy functions and energy density as the 
expansion-normalized bulk viscosity coefficient is increased while the expansion-normalized shear viscosity coefficient is assumed 
to be negligible. One sees that not only in each case is the energy density bounded such that > 0, but also £± — ► 0, so the 
model does indeed isotropize. 
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VIII. DISCUSSION OF RESULTS 

First looking at Figs. ([VII Ap and (|VIIB[) . one can see that the cosmological model does not isotropize in all 
directions. In addition, the asymptotic limit of anisotropy in terms of the "freczing-in" of the dynamical variables was 
also noted by Hervik, van den Hoogen and Coley when studying the Bianchi type IV tilted perfect- fluid model [f|. 
It is interesting to note that as the shear viscosity coefficient r/o increased, the anisotropy in the cosmological model 
increased as well. 

However, as can be seen from looking at Figs. (jVII C|) and (|VII Dj) where the bulk viscosity coefficient £o was 
significantly larger than the shear viscosity coefficient, the model isotropized. What is additionally interesting is that 
as the bulk viscosity was increased, the relationship between the dynamical variables became increasingly linear, which 
hints at the notion that the dynamics of the transitional period from the early-universe cosmology to the isotropic 
universe models become simpler. This provides evidence that the presence of a significant bulk viscous pressure 
rather than a significant shear viscosity is necessary for such a cosmological model to isotropize. These notions were 
also confirmed in the plots of the limiting cases in Figs. (jVII E[) and (jVII Ep . In Fig. (|VII E[) , the different plots in 
this figure demonstrated the behaviour of the anisotropy functions and energy density as the expansion-normalized 
bulk viscosity coefficient was increased while the expansion-normalized shear viscosity coefficient was assumed to be 
negligible. One saw that not only in each case was the energy density bounded such that il > 0, but also E± — > 0, so 
the model did indeed isotropize. On the other hand, in Fig. (jVII E[) . the different plots demonstrated the behaviour 
of the anisotropy functions and energy density as the expansion-normalized shear viscosity coefficient was increased 
while the expansion-normalized bulk-viscosity coefficient was assumed to be negligible. We observed that although 
in each case the energy density was bounded such that f2 > 0, £± -ft 0, so the model did not isotropize. In fact, we 
additionally observed that the anisotropic variables "freeze in" to non-zero asymptotic values after a period of time. 

Similar results for other Bianchi models have been reported by van den Hoogen and Coley [3] , Singh and Kale [l8| , 
and Pradhan, Rai and Singh (Bj]. It is also important to note that no cosmological model that is spatially homogeneous 
and anisotropic and necessarily has a three-dimensional isometry group can evolve to a model with a six-dimensional 
isometry group. When one speaks of a model isotropizing, £± — > only asymptotically. However, any Bianchi model 
can become arbitrary close to that of the FLKW models as a consequence of asymptotic isotropization. 

IX. CONCLUSIONS 

In this paper, we were interested in formulating a viscous model of the early universe based on The Bianchi Type 
IV algebra. We first began by considering a congruence of fluid lines in spacetime, upon which, analyzing their 
propagation behaviour, we rigorously derived the famous Raychaudhuri equation, suitably modified to accommodate 
viscous fluid effects. We then went through in detail the topological and algebraic structure of a Bianchi Type IV 
algebra, by which we derived the corresponding structure and constraint equations. From this, we looked at the 
Einstein field equations in the context of orthonormal frames, and derived and discussed the resulting dynamical 
equations: The Raychaudhuri Equation, generalized Friedmann equation, shear propagation equations, and a set of 
non-trivial constraint equations. It is shown that for cases where the expansion-normalized bulk viscosity coefficient 
£o dominates the expansion- normalized shear viscosity coefficient, rjo, the model isotropizes. 
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